/*
	Интерполяция сплайнами
*/

#include <math.h>
#include <stdlib.h>
#include <assert.h>
#include "sinpol.h"

#define sinpol_alloc_mem(mem_len) (calloc((mem_len), sizeof(double)))

/* Создание объекта полинома */
sinpol_t *sinpol_new(const int sinpol_degr)
{
	assert(sinpol_degr > 1);

	/* Выделение памяти под полином */
	sinpol_t *sinpol_n = malloc(SINPOL_T_SZ);
	if(!sinpol_n) return 0;

	/* Выделение памяти под данные полинома */
	void *mem_n = 0;

	/* Выделение памяти под значения X */
	mem_n = sinpol_alloc_mem(sinpol_degr);
	if(!mem_n) return 0;

	sinpol_n->vals_x = mem_n;

	/* Выделение памяти под значения Y */
	mem_n = sinpol_alloc_mem(sinpol_degr);
	if(!mem_n) return 0;

	sinpol_n->vals_y = mem_n;

	/* Выделение памяти под коэффициенты сплайна */
	mem_n = calloc(sinpol_degr, SINPOL_TUP_T_SZ);
	if(!mem_n) return 0;

	sinpol_n->splines = mem_n;

	/* Установка параметров полинома */
	sinpol_n->sinpol_degr = sinpol_degr;

	return sinpol_n;
}

/* Удаление объекта полинома */
void sinpol_destroy(sinpol_t *sinpol_o)
{
	assert(sinpol_o);

	free(sinpol_o->vals_x);
	free(sinpol_o->vals_y);
	free(sinpol_o->splines);
	free(sinpol_o);
}

/* Генерирование исходных данных */
void sinpol_gen_data(sinpol_t *sinpol_o, double arg_x_start, double arg_x_end, double (*gen_func) (const double arg_x))
{
	assert(sinpol_o);
	assert(gen_func);

	/* Количество исходных точек */
	int pnts_c = sinpol_o->sinpol_degr;

	/* Интервал значений */
	double arg_x_step = (arg_x_end - arg_x_start) / (double) (pnts_c - 1);

	/* Генерирование данных */
	for(int i = 0; i < pnts_c; i++)
	{
		*((sinpol_o->vals_x) + i) = arg_x_start;
		*((sinpol_o->vals_y) + i) = gen_func(arg_x_start);
		
		arg_x_start += arg_x_step;
	}
}

/* Построение сплайна */
void sinpol_build(sinpol_t *sinpol_o)
{
	assert(sinpol_o);

	/* Количество точек */
	int pnts_c = sinpol_o->sinpol_degr;

	/* Инициализация массива сплайнов */
	for(int i = 0; i < pnts_c; ++i)
	{
		sinpol_o->splines[i].arg_x = sinpol_o->vals_x[i];
		sinpol_o->splines[i].arg_a = sinpol_o->vals_y[i];
	}
	sinpol_o->splines[0].arg_c = sinpol_o->splines[pnts_c - 1].arg_c = 0;

	/* Выделение памяти под прогоночные коэффициенты */
	double *args_alpha = sinpol_alloc_mem(pnts_c - 1);
	assert(args_alpha);

	double *args_beta = sinpol_alloc_mem(pnts_c - 1);
	assert(args_beta);

	args_alpha[0] = args_beta[0] = 0;

	/* Решение СЛАУ */
	double *vals_x = sinpol_o->vals_x;
	double *vals_y = sinpol_o->vals_y;

	for(int i = 1; i < pnts_c - 1; ++i)
	{
		double hi_f = vals_x[i] - vals_x[i - 1];
		double hi_s = vals_x[i + 1] - vals_x[i];
		double arg_a = hi_f;
		double arg_b = hi_s;
		double arg_c = 2 * (hi_f + hi_s);
		double arg_f = 6 * ((vals_y[i + 1] - vals_y[i]) / hi_s - (vals_y[i] - vals_y[i - 1]) / hi_f);
		double arg_z = arg_a * args_alpha[i - 1] + arg_c;

		args_alpha[i] = -arg_b / arg_z;
		args_beta[i] = (arg_f - arg_a * args_beta[i - 1]) / arg_z;
	}

	/* Поиск решения (обратный ход метода прогонки) */
	for (int i = pnts_c - 2; i > 0; --i)
	{
		sinpol_o->splines[i].arg_c = args_alpha[i] * sinpol_o->splines[i + 1].arg_c + args_beta[i];
	}

	/* Поиск коэффициентов B и D */
	for(int i = pnts_c - 1; i > 0; --i)
	{
		double hi_f = vals_x[i] - vals_x[i - 1];
		
		sinpol_o->splines[i].arg_d = (sinpol_o->splines[i].arg_c - sinpol_o->splines[i - 1].arg_c) / hi_f;

		sinpol_o->splines[i].arg_b = hi_f * (2 * sinpol_o->splines[i].arg_c + sinpol_o->splines[i + 1].arg_c) / 6 + (vals_y[i] - vals_y[i - 1]) / hi_f;
	}

	/* Освобождение памяти */
	free(args_alpha);
	free(args_beta);
}

/* Вычисление значения в точке */
double sinpol_get_val(const sinpol_t *sinpol_o, const double arg_x)
{
	assert(sinpol_o);

	int pnts_c = sinpol_o->sinpol_degr - 1;
	sinpol_tup_t sinpol_tup_o;

	if(arg_x <= sinpol_o->splines[0].arg_x)
	{
		sinpol_tup_o = sinpol_o->splines[0];
	}
	else if(arg_x >= sinpol_o->splines[pnts_c - 1].arg_x)
	{
		sinpol_tup_o = sinpol_o->splines[pnts_c - 1];
	}
	else
	{
		int i = 0;
		int j = pnts_c - 1;

		while(i + 1 < j)
		{
			int k = i + (j - i) / 2;

			if(arg_x <= sinpol_o->splines[k].arg_x)
			{
				j = k;
			}
			else
			{
				i = k;
			}
		}

		sinpol_tup_o = sinpol_o->splines[j];
	}

	double dx = arg_x - sinpol_tup_o.arg_x;
	double res = sinpol_tup_o.arg_a + (sinpol_tup_o.arg_b + (sinpol_tup_o.arg_c / 2 + sinpol_tup_o.arg_d * dx / 6) * dx) * dx;

	return res;
}